################################################################################
# Function to Calculate ACME  

ACME.fun <- function(m.comp, y.comp){
  
  ACME.mat <- matrix(NA, 10000, 3)
  
  X.m <- m.comp$model
  X.m[,1] <- rep(1, dim(X.m)[1])
  X0.m <- X.m
  X1.m <- X.m
  X0.m[,2] <-  rep(0, dim(X.m)[1])
  X1.m[,2] <-  rep(1, dim(X.m)[1])
  
  for(s in 1:10000){
    ### Sim 
    m.sim <- sim(m.comp, 1)
    
    b.m <- coef(m.sim)[1,]
    m.hat.0 <- b.m%*%t(X0.m)  
    m.hat.1 <- b.m%*%t(X1.m)
    
    X.y <- y.comp$model
    X.y[,1] <- rep(1, dim(X.y)[1])
    X0.y <- X.y
    X1.y <- X.y
    X0.y[,2] <-  rep(0, dim(X.y)[1])
    X1.y[,2] <-  rep(1, dim(X.y)[1])
    
    X00.y <- X0.y
    X00.y[,3] <- t(m.hat.0)
    X00.y <- cbind(X00.y, t(m.hat.0)*0)
    
    X01.y <- X0.y
    X01.y[,3] <- t(m.hat.1)
    X01.y <- cbind(X01.y, t(m.hat.1)*0)
    
    X10.y <- X1.y
    X10.y[,3] <- t(m.hat.0)
    X10.y <- cbind(X10.y, t(m.hat.0)*1)
    
    X11.y <- X1.y
    X11.y[,3] <- t(m.hat.1)
    X11.y <- cbind(X11.y, t(m.hat.1)*1)
    
    y.sim <- sim(y.comp, 1)
    b.y <- coef(y.sim)[1,]
    
    y.hat.00 <- b.y%*%t(X00.y)
    y.hat.01 <- b.y%*%t(X01.y)
    y.hat.10 <- b.y%*%t(X10.y)
    y.hat.11 <- b.y%*%t(X11.y)
    
    ACME.mat[s,1] <- mean(y.hat.01 - y.hat.00)
    ACME.mat[s,2] <- mean(y.hat.11 - y.hat.10)
    ACME.mat[s,3] <- mean(((y.hat.01 - y.hat.00)+(y.hat.11 - y.hat.10))/2)
  }
  
  return(ACME.mat)
}

